% compute f(x) = max_i w_i Gauss(x | mu(i), Sigma(i))
% Code based on 
%http://staff.itee.uq.edu.au/marcusg/msg.html

function maxGMMplot()

rng(42) % seed
params = initialize (2, 20, 5, -5, 10, 0.8);
plotlandscape (5, -5, 50, params); 

end

function params =  initialize(dimension,nGaussian,upper,lower,globalvalue,ratio)

%--------------------------------------------------------------------------
%This is the initialization function of the Gaussian landscape generator
%
%Syntax: initialize(dimension,nGaussian,upper,lower,globalvalue,ratio)
%
%Example:initialize(2,10,5,-5,1,0.8)
%
%Inputs: 
%        dimensionality
%        number of Gaussian components
%        upper and lower boundaries
%        value of the global optimum
%        values of local optima ([0,ratio*globalvalue])
%
%Outputs:
%        inverse covariance matrix
%        meanvector
%        component peak value
%
%Author: Bo Yuan (boyuan@itee.uq.edu.au)
%--------------------------------------------------------------------------

%{
clear global covmatrix_inv;
clear global meanvector;
clear global optimumvalue;

global covmatrix_inv;   %the inverse covariance matrix of each component
global meanvector;      %the mean of each component
global optimumvalue;    %the peak value of each component
%}



if nargin<6
    
    disp('Usage: initialize(dimension,nGaussian,upper,lower,globalvalue,ratio)');
    return;    
    
end

if dimension<=1|nGaussian<=0|upper<=lower|globalvalue<=0|ratio<=0|ratio>=1
    
    disp('Incorrect parameter values!');
    return;
    
end

% Generate rotation matrix

e=diag(ones(1,dimension));   % unit diagonal matrix

for i=1:nGaussian
    
   rotation{i}=e;            % initial rotation matrix for each Gaussian
   
end

for i=1:nGaussian            

  for j=1:dimension-1        % totally n(n-1)/2 rotation matrice
      
      for k=j+1:dimension
          
        r=e;
        
        alpha=rand*pi/2-pi/4;% random rotation angle [-pi/4,pi/4]
        
        r(j,j)=cos(alpha);
        r(j,k)=sin(alpha);
        r(k,j)=-sin(alpha);
        r(k,k)=cos(alpha);
    
        rotation{i}=rotation{i}*r;
    
     end
    
  end
  
end

% Generate covariance matrix

variancerange=(upper-lower)/20;   % this controls the range of variances

variance=rand(nGaussian,dimension)*variancerange+0.05;  % add 0.05 to avoid zero variance

for i=1:nGaussian

      covmatrix=diag(variance(i,:));
      covmatrix=rotation{i}'*covmatrix*rotation{i};
      covmatrix_inv{i}=inv(covmatrix);
      
end

% Generate mean vectors randomly within [lower, upper]
meanvector=rand(nGaussian,dimension)*(upper-lower)+lower; 

% assign values to components
optimumvalue(1)=globalvalue;     % the first Gaussian is set to be the global optimum

% values of others are randomly generated within [0,globalvalue*ratio]
optimumvalue(2:nGaussian)=rand(1,nGaussian-1)*globalvalue*ratio;

params.covmatrix_inv = covmatrix_inv;
params.meanvector = meanvector;
params.optimumvalue = optimumvalue;

end


function [fitnessvalue,components]=fitness(x, params)

%--------------------------------------------------------------------------
%This is the fitness function of the Gaussian landscape generator
%
%Syntax: [fitnessvalue,components]=fitness(x)
%
%Inputs:
%        individuals stored in a p-by-n matrix x
%
%Outputs:
%        the fitness value of each row in x
%        the value of each row in x generated by each component
%
%Author: Bo Yuan (boyuan@itee.uq.edu.au)
%--------------------------------------------------------------------------

%{
global covmatrix_inv;  %the inverse covariance matrix of each component
global meanvector;     %the mean of each component
global optimumvalue;   %the peak value of each component

if isempty(covmatrix_inv)|isempty(meanvector)|isempty(optimumvalue)
    
    disp('Run initialize function first!');
    return;
    
end
%}


covmatrix_inv = params.covmatrix_inv;
meanvector = params.meanvector;
optimumvalue = params.optimumvalue;


nGaussian=size(meanvector,1);  % total number of components

[p,n]=size(x);                 % p: number of individuals; n: dimensionality

if n~=size(meanvector,2)
    
    disp('Incorrect dimension of inputs!');
    return;
    
end

tmp=zeros(nGaussian,p);

%----------------------------------------------------

for i=1:nGaussian              % calculate the values generated by each component
    
    newx=x-repmat(meanvector(i,:),p,1);
    
    z=(newx*covmatrix_inv{i}).*newx;
           
    tmp(i,:)=sum(z,2)';        
   
end

f=exp(-0.5*tmp/n);             % f is a nGaussian-by-p matrix

f=f.*repmat(optimumvalue',1,p);% multiply the peak value of each component

components=f';                 % the value of each individual generated by each component

fitnessvalue=max(f,[],1);      % choose the maximum values as the fitness values
end



function plotlandscape(upper,lower,N, params)

%--------------------------------------------------------------------------
%This is the plotting function of the Gaussian landscape generator(2D only)
%
%Syntax: plotlandscape(upper,lower,N)
%
%Example:plotlandscape(5,-5,100)
%
%Inputs:
%       upper boundary of the search space
%       lower boundary of the search space
%       number of sampling points
%
%Outputs:
%       3D surface plot
%       2D contour plot
%
%Author: Bo Yuan (boyuan@itee.uq.edu.au)
%--------------------------------------------------------------------------


if upper<=lower
    
    disp('Upper must be larger than Lower!');
    return;
    
end

if N<10
    
    disp('Incorrect N value!');
    return;
    
end

%-----------------------------------------------------

inc=(upper-lower)/N;

x=lower:inc:upper;     %x coordinates
y=lower:inc:upper;     %y coordinates

d=length(y);

pos=zeros(d*d,2);      %(x,y)coordinates for all sampling points

for i=1:d
    
    pos((i-1)*d+1:i*d,1)=x(i)*ones(d,1);
    
end

pos(:,2)=repmat(y',d,1);

f=fitness(pos, params);        %evaluate individuals

z=vec2mat(f,d)';       %transform into a matrix
figure;
  
colormap(gray);        %3D surface plot
surfl(x,y,z);
shading interp;

figure;

[C,H]=contour(x,y,z);  %2D contour plot
end

% https://searchcode.com/codesearch/view/9557982/
function [M, d] = vec2mat(V, c, val)

  switch (nargin)
    case 1,
      M = V;
      return;
    case 2,
      val = 0;
    case 3,
      val = val;
    otherwise
      error ("usage: [M, add] = vec2mat (V, c [, d])");
  end

  V = V.';
  V = V(:);

  r = ceil (length (V) / c);

  d = r * c - length (V);
  if (d ~= 0)
    V = [ V ; val*ones(d, 1) ];
  end

  M = reshape (V, c, r).';

  end

